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ABSTRACT 

A  fully  Sinc-Galerkin  method  in  both  space  and  time  is  presented  for  fourth-order  time- 
dependent  partial  differential  equations  with  fixed  and  cantilever  boundary  conditions.  The 
sine  discretizations  for  the  second-order  temporal  problem  and  the  fourth-order  spatial  prob¬ 
lems  are  presented.  Alternate  formulations  for  variable  parameter  fourth-order  problems  are 
given  which  prove  to  be  especially  useful  when  applying  the  forward  techniques  of  this 
paper  to  parameter  recovery  problems.  The  discrete  system  which  corresponds  to  the  time- 
dependent  partial  differential  equations  of  interest  are  then  formulated.  Computational  issues 
are  discussed  and  a  robust  and  efficient  algorithm  for  solving  the  resulting  matrix  system  is 
outlined.  Numerical  results  which  highlight  the  method  are  given  for  problems  with  both 
analytic  and  singular  solutions  as  well  as  fixed  and  cantilever  boundary  conditions. 
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1  Introduction 


The  Sinc-Galerkin  method  for  partial  differential  equations  (PDE’s)  has  previously  been 
developed  for  the  model  elliptic  problem  in  two  and  three  dimensions  [1,  2],  the  parabolic 
problem  in  one  and  two  dimensions  [3,  2],  and  the  second-order  hyperbolic  problem  in  one 
dimension  [4].  The  present  work  extends  the  method  to  fourth-order  time- dependent  prob¬ 
lems  with  various  common  boundary  conditions.  This  extension  is  important  for  the  very 
practical  reason  that  the  numerical  solution  of  problems  in  this  class  is  necessary  in  ap¬ 
plications  ranging  from  the  control  of  large  flexible  space  structures  to  the  development  of 
robotics  designs  [5,  6,  7]. 

For  clarity  of  development,  the  method  will  be  presented  for  the  linear  fourth-order  time- 
dependent  problems 

d^xi  d ^  (  3^u  \ 

Cu(x,t)  =  -QP  (x,t)  +  —  (  EI(x)-^(x,t)J  =  f(x,t),  0<®<1  *>0 

u(0,t)  =  u(l,t)  =  0,  t>  0  (1.1) 

|(« .0- £(».') -o,  ‘>o 

du 

u(i,0)  =  -^(x,0)  -  0,  0  <  x  <  1 

and 

Cu(x,t)  =  f{x,t),  0  <  x  <  1  t  >  0 

“(0,0  =  5(0,  (1,0  =  7(0.  <>0  (1.2) 

^(o, 0  =  3(0,  s(«^)d.0-I( 0,  <>o 

du 

u(x,0)  =  —(x,0)  =  0,  0  <  x  <  1. 

These  formulations  are  generalizations  of  the  equations  which  arise  when  using  the  Euler- 
Bernoulli  theory  to  model  beams  with  flexural  rigidity  EI( x)  and  fixed  and  cantilevered 
ends,  respectively.  The  general  j(t)  and  6(t )  in  (1.2)  allow  for  the  inclusion  of  boundary 
controllers.  For  ease  of  presentation,  the  boundary  conditions  in  (1.1)  and  (1.2)  will  respec¬ 
tively  be  referred  to  as  fixed  and  cantilever  conditions  throughout  the  paper.  It  is  noted  that 
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the  methods  of  this  work  are  easily  extended  to  problems  with  simple  and  free  boundary 
conditions  with  further  details  given  in  [8]. 

The  construction  of  an  approximate  solution  to  problems  of  the  form  (1.1)  or  (1.2)  com¬ 
monly  begins  with  a  Galerkin  discretization  of  the  spatial  variable  with  time-dependent 
coefficients.  This  yields  a  system  of  ordinary  differential  equations  which  is  solved  via  dif¬ 
ferencing  techinques.  Due  to  stability  constraints  on  the  discrete  evolution  operator,  low 
order  methods  with  small  time  steps  are  often  required  to  obtain  accurate  approximations. 
In  contrast,  the  method  of  this  work  implements  a  Galerkin  scheme  in  time  as  well  as  space. 
Because  the  basis  functions  are  tensor  products  of  sine  functions  composed  with  suitable 
conformal  maps,  the  method  has  the  inherent  advantage  that  the  study  of  error  analysis  and 
matrix  structure  begins  at  the  level  of  an  ordinary  differential  equation. 

The  fully  Sinc-Galerkin  method  in  space  and  time  has  many  other  salient  features  due 
both  to  the  properties  of  the  basis  functions  and  to  the  manner  in  which  the  problem  is 
discretized.  First,  the  judicious  choice  of  a  conformal  map  provides  approximate  solutions 
to  (1.1)  and  (1.2)  which  are  valid  on  the  infinite  time  interval  rather  than  only  on  a  trun¬ 
cated  time  domain.  Furthermore,  the  optimal  exponential  convergence  rate  is  maintained 
even  in  the  presence  of  boundary  singularities.  Finally,  the  discrete  system  requires  no  nu¬ 
merical  integrations  to  fill  either  the  coefficient  matrix  or  the  right-hand  side  vector.  All 
three  features  prove  to  be  advantageous  when  solving  both  forward  and  inverse  fourth-order 
problems.  A  drawback  to  the  method  is  that  it  produces  a  full  system  in  contrast  to  the 
banded  matrices  which  are  associated  with  finite  difference  and  finite  element  methods.  In 
part,  the  exponential  convergence  rate  offsets  this  disadvantage. 

The  foundations  of  the  Sinc-Galerkin  method  are  described  in  Section  2.  The  fundamental 
quadrature  rule  is  given,  and  the  exponential  convergence  rate  of  this  method  is  stated.  A 
thorough  review  of  sine  function  properties  can  be  found  in  [9]  and  [10]. 

In  the  next  section,  the  discretization  of  the  second-order  temporal  and  fourth-order  spa¬ 
tial  operators  is  outlined.  Two  schemes  for  discretizing  EI(x)  are  presented.  These  two 
schemes  are  motivated  on  the  one  hand  by  the  forward  problem  and  on  the  other  hand  by 
the  inverse  problem  involving  the  recovery  of  EI(x),  given  sampled  data.  In  the  first,  EI(x) 
is  differentiated  directly,  whereas  in  the  second  (as  motivated  by  the  parameter  recovery 
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problem)  EI(x)  is  replaced  by  a  finite  dimensional  expansion  EImm  before  differentiation. 
Attention  is  focussed  on  preserving  the  method’s  exponential  convergence  rate  while  dis¬ 
cretizing  EI(x )  and  adapting  to  varying  boundary  conditions. 

In  Section  4,  the  one-dimensional  results  from  the  previous  section  are  combined  to  yield 
methods  which  very  accurately  approximate  the  solutions  to  the  fourth-order  time-dependent 
problems  (1.1)  and  (1.2).  Two  equivalent  formulations  for  the  resulting  discrete  system  are 
presented  and  a  very  robust  and  accurate  solution  algorithm  is  outlined. 

Numerical  results  are  presented  in  the  fifth  section.  Of  the  many  examples  tested,  those 
difcussed  in  this  section  best  exhibit  the  features  necessary  for  the  practical  implementation 
of  the  Sine- Galerkin  method.  The  first  and  second  examples  illustrate  the  method  as  applied 
to  problems  with  fixed  boundary  conditons  while  the  third  and  fourth  examples  have  can¬ 
tilever  boundary  conditions.  The  first  and  last  examples  have  analytic  solutions,  the  second 
example  has  an  algebraic  singularity,  and  the  third  example  contains  a  logarithmic  singular¬ 
ity.  The  numerical  results  demonstrate  that  the  exponential  convergence  rate  is  maintained 
in  all  four  cases. 
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□  □ 


Figure  1.  Three  Adjacent  Members  (S(k,  h*)(x),  k  =  —1,0,1,  h*  =  of  the  Translated  Sine 
Family. 

To  construct  basis  functions  on  the  intervals  (0,1)  and  (0,oo),  respectively,  consider  the 
conformal  maps 

m  =  In  (j-^)  (2.1) 

and 

T(u/)  =  ln(w).  (2.2) 

The  map  <f>  carries  the  eye- shaped  region 

De  -  =  *  +  iy  :  j arg  ( |  <  d  < 

onto  the  infinite  strip 

Ds  =  =  C+iv  •  \v\  <  d  <  ^}- 

Similarly,  the  map  T  carries  the  infinite  wedge 

Dw  —  jttf  =  t  +  is  :  |ar(7(u;)|  <  d  <  —  j 

onto  the  strip  D$.  These  regions  are  depicted  in  Figure  2. 


4 


Figure  2.  The  Domains  Ds,  De,  and  Dw- 


The  sine  gridpoints  z k  €  (0,1)  in  De  will  be  denoted  x*  since  they  are  real.  Similarly,  the 
gridpoints  i vk  £  (0,oo)  in  Dw  will  be  denoted  tk.  Both  are  inverse  images  of  the  equispaced 
grid  in  Ds\  that  is 


Xfc  =  <f>  i(kh)  = 


e 


kh 


1  +  ekh 


and 


tk  =  T~l(kh)  =  ekh. 


To  simplify  notation  throughout  the  remainder  of  this  section,  the  pairs  De  and  T,  Dw 
are  referred  to  generically  as  x>  D.  It  is  understood  that  the  subsequent  definition  and 
theorems  hold  in  either  setting.  Furthermore,  the  inverse  of  x  is  denoted  by  V’- 

The  important  class  of  functions  for  sine  interpolation  and  quadrature  is  denoted  B(D ) 
and  defined  next. 
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Definition  2.1.  Let  B(D )  be  the  class  of  functions  F  which  are  analytic  in  D,  satisfy 


f  \F(z)dz\  — >  0,  t  —*  ioo 

where  L  =  {is  :  js]  <  d  <  *},  and  on  the  boundary  of  D  (denoted  dD)  satisfy 

N{F)  =  f  \F(z)dz\  <  oo. 

J8D 

The  following  theorem  for  functions  in  B(D)  is  found  in  [11]. 


Theorem  2.1.  Let  T  be  (0,1)  or  (0,oo)  when  x  —  <t>  or  T,  respectively.  If  F  £  B(D)  and 
Zj  =  ip(jh)  =  =  0,±1,±2,...,  then  for  h  >  0  sufficiently  small 


fF(z)dz-h  £ 


n*i) 


]=-00 


X'(*j) 


<  Kxe 


-lird/h 


(2.3) 


Theorem  2.1  illustrates  the  exponential  convergence  rate  which  is  a  trademark  of  sine  meth¬ 
ods.  There  is  a  common  occasion  when  it  is  possible  to  evaluate  the  infinite  series  appearing 
in  (2.3),  namely  when  integrating  against  S(k,  h)  o  x ■  In  general,  however,  the  series  must 
be  truncated.  With  additional  hypotheses,  it  is  proven  in  [9]  and  [12]  that  the  truncation 
need  not  be  at  the  expense  of  the  exponential  convergence. 


Theorem  2.2.  Assume  F  £  B(D )  and  that  there  exist  positive  constants  K,a  and  (3  such 
that 

( 

exP(-Q!x(T)l)>  t  €  t/>((-oo,0)) 
exP(— /2|x(t)I)>  t£^([0,oo)). 


F(t) 


|x'(t) 

Then  for  h  sufficiently  small 


<  K 


(2  4! 

V4-' 


<  Kxe~2lrJ^h  +  - e-aMh  +  -e~0Nh. 
a  13 


(2.5) 


Theorems  2.1  and  2.2  are  used  to  establish  a  uniform  error  bound  when  building  an 
approximate  solution  to  an  ordinary  differential  equation  (ODE).  It  should  be  noted  that 
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the  nature  of  the  class  B(D)  guarantees  that  the  exponential  convergence  rate  holds  for 
many  differential  equations  with  singular  solutions;  that  is,  problems  where  the  solution  has 
an  unbounded  derivative  on  the  boundary.  By  applying  the  scheme  to  select  second-  and 
fourth-order  ODE’s,  one  can  derive  the  fundamental  matrices  comprising  the  discrete  sine 
system  for  the  fourth-order  time-dependent  problems  of  interest. 


3  Sinc-Galerkin  Systems  for  ODE’s 


In  this  section,  the  sine  discretizations  will  be  catalogued  for  the  second-order  temporal 
problem  and  two  different  fourth-order  spatial  problems  distinguished  by  their  boundary 
conditions.  Alternate  formulations  for  the  variable-parameter  fourth-order  problems  will 
also  be  given  which  prove  to  be  especially  useful  when  applying  the  forward  techniques 
outlined  in  this  paper  to  parameter  recovery  problems. 

In  order  to  construct  the  discrete  Sinc-Galerkin  system  for  either  the  temporal  or  spatial 
problems,  the  following  identities  are  needed.  Let 


2  [S(p,  h)  o  x(2)l  = 

Z-Zi 


1,  i  =  p 

0,  i^P, 


^]  =  h  -j-S{p,h)oX{z) 

■  X  t=z{ 

4?  -  h*  -fis(p>h)°x(z) 

■  x  z-z. 


(-I)*-" 

(*-p) 


l  =  p 


,  *  ^  p, 


3  ’ 

(— 2)(-l)‘-p 


z  =  p 


(*  ~P) 


—  > 


$  =  ^  JlS(P>h)°x{z) 

■  X  z=Zi 


0, 

(-iy-p 

(f  -  P)3 


* =p 


[6  —  7TJ(t  —  p)2],  i^p, 
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and 


f  tt4 

1  =  P 

*„(_!)<->>  (3-5) 

denote  the  evaluation  at  the  gridpoint  z,  of  the  sinc-map  compositions  and  their  derivatives 
with  respect  to  the  map  x- 

S.l.  The  Spatial  Problem:  Fixed  Boundary  Conditions 

In  [13],  a  thorough  analysis  of  the  Sinc-Galerkin  method  is  given  for  linear  fourth-order 
ODE’s  with  fixed  boundary  conditions.  For  purposes  of  constructing  the  sine  discretization 
for  (1.1),  it  suffices  to  review  that  procedure  for 


-  * 


dA 


—S{p,h)oX(z) 


Lu(x)  =  ( EI(x)u"(x ))"  =  fix),  0  <  x  <  1 

u(0)  =  it(l)  =  0  (3-6) 

u'(0)  =  ti'(l)  =  0. 

Note  that  the  interval  (0,1)  is  for  convenience  only;  adapting  the  map  <f>  (see  (2.1))  generalizes 
the  method  to  any  finite  interval  (o,  b). 

To  define  the  Sinc-Galerkin  approximation  to  (3.6),  select  the  basis  {5’i}^_Ma.  where 
Si{x)  =  5(r,  hx)  o  <f){x)  and  take  the  approximate  solution  to  be 

AT. 

«m.(*)  =  U»5'*(X)>  mx  -  Mx  +  Nm  +  1.  (3.7) 

i=-M. 

The  unknown  coefficients  {u*}  in  (3.7)  are  determined  by  orthogonalizing  the  residual 
lum,  -  /  with  respect  to  the  functions  This  yields  the  discrete  system 

iLumm  —  /,  Sp)  —  0 


for  p  =  —  Mx,  •  •  • ,  Nt.  The  weighted  inner  product  (•,  ■)  is  taken  to  be 

(F,  G)  =  f  F{x)G(x)w{x)  dx 
Jo 


where 


w(x)  = 


1 

(Si? 


(3.8) 

(3.9) 
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For  a  further  discussion  concerning  the  choice  of  weight,  see  [13]. 

Before  invoking  the  quadrature  rules,  integration  by  parts  is  used  to  transfer  the  differ¬ 
entiation  of  u  onto  Spw ,  thus  yielding  the  system 

/  tt(x)[.EJ(x)(Sp(x)iu(x))'/]/'dx  =  /  f(x)Sp(x)w(x)dx  (3.10) 

Jo  Jo 

for  p  =  —  MX)  •  •  • ,  Nx,  With  the  weight  choice  (3.9),  the  boundary  terms 

{(EIu")'{Spw)  -  (EIu")(Spw)'  +  u'(EI(Spw)')'  -  u(EI(Spw)')"}(x) '  (3.11) 

o 

vanish  for  essentially  all  problems  of  interest. 

Two  approaches  are  distinguished  by  the  treatment  of  the  first  integral  in  (3.10).  In  the 
traditional  scheme,  EI(x)  is  differentiated  directly  and  the  resulting  integrals  are  approx¬ 
imated  via  sine  quadrature  rules.  This  scheme  is  direct  and  suitable  for  a  large  class  of 
forward  problems.  The  second,  alternative,  approach  is  motivated  by  the  parameter  recov¬ 
ery  problem  and  differs  from  the  first  in  that  EI(x)  is  replaced  by  a  sine  expansion,  EImm, 
before  quadrature  is  applied.  Both  approaches  then  proceed  in  the  same  manner  whereby 
the  system  is  expanded  and  the  resulting  integrals  are  evaluated  via  Theorem  2.2,  or  when 
possible,  Theorem  2.1. 

The  careful  choice  of  the  decay  parameters  a  and  (3  in  (2.4)  provides  a  means  of  balancing 
the  asymptotic  errors  resulting  from  the  quadrature  and  hence  minimizes  the  system  size. 
With  regard  to  (2.4),  the  condition 

|E/(i)ti(i)u>(i)(#'i (x))3|  <  K  (  2)  (3.12) 

guarantees  the  decay  needed  to  truncate  the  infinite  quadrature  rule.  A  less  general  but 
more  convenient  assumption  than  (3.12)  is 

j£/(x)tt(x)j  <  Kxa+*{  1  -  x)^*.  (3.13) 


With  a  and  /?  specified  and  M*  chosen,  the  parameter  selections 


(3.14) 
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and 


~MX  +  1 

/3 


(3.15) 


balance  the  asymptotic  quadrature  errors  to  at  least  ).  This  rate  results  from 

the  presence  of  a  sine  function  in  the  integral.  In  the  above,  [•]  denotes  the  greatest  integer 
function.  Note  that  if  jMx  is  an  integer,  (3.15)  can  be  replaced  by  the  selection  Nx  =  %MX. 

The  discrete  system  for  (3.6),  using  the  traditional  approach,  can  then  be  formulated 
as  follows.  Let  I^l\i  =  0, 1,2, 3, 4  denote  the  mx  x  mx  matrices  whose  pi-th  entry  is 
from  (3.1)  -  (3.5)  and  let  V{t})  be  the  diagonal  matrix  with  entries  7?(i_m,,)>  •  •  • ,  v{xNm)- 
The  vector  of  unknowns  u  =  •  •  • , «n„]T  is  then  related  to  the  known  vector  /  = 


[/(z-aO.  —  ./Oew.)]7'  by 


AxH=V{{<j>')-')f 


(3.16) 


where 


-^/(4>r>(a4)  4-  ~I(3)V(a3)  +  ^IWV(a2)  +  0  +  I^V(a0) 

tix  hzx  .  hx 


The  functions  at(x),l  =  0, 1,2, 3, 4  are  given  by 


ao 


(. EIw )""  „( El'w )'"  ( EI"w )" 

4>>  4>'  +  4 


AH  AW  AWt 

ax  =  A(EIw)'"  +  G(EIw)"~  +  4(£/u>)'^-  +  EIw?— 

t pf  &  (T)' 

AH  lW  AH 

-6{EI'w)"  -  6 {EI'w)'^-  -  2 El'w?-  +  El"{w)^~  +  2 {EI"w)\ 

<P  <P  <t> 

a2  =  6(EIw)"0'  +  12  {EIw)'4>"  +  4  EIw<f>'"  +  3  (EIw)'^^- 

<t>' 

-6(E/'tx;)V'  -  SEIw<j>"  +  El"w4>\ 

a3  =  4 {EIw)\4>'Y  +  6EIw<f>'<t>"  -  2 El'w{4>'f 


and 


o4  =  EIw(4>')\ 


(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 


Further  details  concerning  the  derivation  of  the  system  (3.16)  and  a  thorough  spectral  anal¬ 
ysis  of  the  component  matrices  can  be  found  in  [13]. 
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As  mentioned  previously,  the  treatment  of  the  first  integral  in  (3.10)  yields  various  pertur¬ 
bations  of  the  method  which  are  advantageous  in  certain  applications.  One  such  application 
is  the  parameter  recovery  problem  where  an  integral  part  of  most  numerical  schemes  for 
solving  that  problem  is  an  accurate  forward  solver.  With  this  in  mind,  the  alternative  ap¬ 
proach  mentioned  above  is  implemented  wherein  the  term  EI(x)  in  (3.10)  is  expanded  as  a 
linear  combination  of  weighted  sine  functions  with  four  Hermite-like  algebraic  terms.  These 
terms  are  added  to  accommodate  the  potentially  nonzero  function  and  derivative  values  of 
El  at  x  =  0  and  x  =  1.  Specifically,  this  parameter  basis  is  taken  to  be  { with 

k  =  -Mx 

b-Mm+ i(a),  fc  =  -M*  +  1 

Mx)='  vE(x)Sk{x),  -Mx  +  2  <  k  <  Nx  -  2  (3-23) 

t>Nm- 1(®),  k  =  Nx-  1 

bNm(x),  k  =  Nx. 

Here  Sk(x)  =  S(k,hx)  o  ^>(x)  and  the  basis  weight  vE  is  taken  to  be 

vE(x)  =  w(x)  =  [x(l  -  x)]$.  (3.24) 

The  algebraic  boundary  basis  functions  are  given  by 

b-M.+i(x)  =  (1  -  x)2[2x  !  1], 
bN.~ l(x)  =  l2[2(l  -  x)  +  1], 

b-M.(x )  =  *0  -  xf 

and 

bN.{x)  =  -x2(l  -  x). 

The  finite  dimensional  approximation  of  El  then  takes  the  form 

£/m.(*)=  ^  CfcV'fc(x).  (3.25) 

k  =  -M. 
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The  number  of  basis  functions  used  in  the  expansion  is  chosen  so  as  to  guarantee  a  square 
coefficient  matrix.  This  is  done  to  simplify  the  implementation  of  the  method  when  applied 
to  the  PDE  (1.1)  of  interest. 

A  quick  note  should  be  made  concerning  the  choice  of  basis  and  the  manner  of  expand¬ 
ing  EIma.  The  two  derivative-interpolating  boundary  basis  functions  are  added  so  that  this 
expansion  of  EImm  is  the  same  as  that  used  with  cantilever  or  free  boundary  conditions.  The 
choice  of  (3.24)  for  basis  weight  is  certainly  sufficient  and  proves  to  be  beneficial  when  in¬ 
corporating  this  forward  scheme  into  a  numerical  method  for  solving  the  parameter  recovery 
problem. 

The  expansion  (3.25)  is  substituted  into  (3.10)  and  the  resulting  integrals  are  evaluated 
via  Theorem  2.2  or  Theorem  2.1  when  possible.  The  decay  condition  (2.4)  equates  to  the 
condition 

* 

£“+5  xe (o  -) 

\Sl{x)u(x)\  <  K  ’  V  ’ 2J  (3.26) 

(1-xf+f,  xe(±,l) 

where  the  “homogeneous”  part  of  El  is 

£J(x)  =  EI(x)  -  EI(0)b.uMx)  -  £/(1)&at.-i(*)  ~  EJ\0)b.Mm{x)  -  EI'{l)bN,(x).  (3.27) 

The  arguments  leading  to  this  condition  are  analogous  to  those  presented  in  the  second-order 
case  as  described  in  [14].  Again,  this  may  be  replaced  by  the  more  stringent  requirement 

|£I(x)u(x)|  <  ffxa+’(l  —  x)0+*. 

As  before,  the  asymptotic  errors  are  balanced  by  choosing  hx  and  Nx  as  specified  in  (3.14) 
and  (3.15). 

With  t2  and  /  defined  as  before  and  El  expanded,  the  system  for  (3.6)  using  this  alter¬ 
native  approach  can  be  written  as 

AxU  =  V({<f,’)-^)f  (3.28) 

where 

Ax  =  [$(J)D(p>))  +  2*<3>P(p„(l))  4-  ^V(p^)).  (3.29) 
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The  notation  T>(p^t))tt  =  0, 1,2  denotes  the  diagonal  matrices  containing  the  components 
of  the  vectors 

=  tf(0c,  £  =  0,1,2  (3.30) 

where  c  =  [c_Mm,  •  •  • ,  cjvJt.  The  matrices  =  2,3,4  and  =  0,1,2  are  defined 

componentwise  by 

[*wW  =  ^(sp®)(<>(*i)  <3-31) 

and 

[*(0U  =  1>P(xi).  (3.32) 


The  notation  on  the  right-hand  sides  of  (3.31)  and  (3.32)  indicates  the  j-th  and  £-th  deriva¬ 
tives,  respectively. 

To  illustrate  the  dependence  of  =  2,3,4  and  'P^,£  =  0, 1,2  on  previously  defined 

matrices,  the  respective  expansions  are  listed  beL^.  The  diagonal  matrices  T>  and  the 
matrices  1^,1  =  0, 1,2, 3, 4  have  sizes  consistent  with  the  following  range  of  the  indices  t,p 
and  k  {—Mx  <  i  <  Ns,  —Mx  <  p  <  Nx,  —Ms  +  2  <  k  <  Nx  —  2).  From  (3.31)  it  follows  that 


and 


$(2)  =  J_/(2)£>(  w<f>')  +  -^-/(1)T>(  2u/)  +  I<0)V(w"), 

“'I 

$(3>  =  +  ~I{2)V{2w'<t>'  +  3 w<f>") 

+t/')v{3w''+3wj+wt)+i1°)v(t) 

$(4)  =  LjWv  [w^'Y)  +  +  Qw<f>'4>") 


(t) 


WT 


+  I2w'<f>"  +  4t u<f>"  +  3u> 


4>' 


iw'"  +  6w"^-  +  4ur-—  +  w~ 


<t>' 


+/(°)p 


For  l  —  0, 1,2  the  mx  x  mx  matrices  in  (3.32)  are  given  by 

*»  -  [a. :  :  b('»  :  <5t. : 


(3.33) 


(3.34) 


(3.35) 


(3.36) 
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where  =  [ij^(x_Af.),  •  •  • ,  &fc^(siO]T  for  k  =  -Mx,  -Mx  +  1,  -N*  —  1,  and  iVx.  Again,  the 
superscript  l  indicates  the  £-th  derivative.  The  mx  x  (mx  -  4)  matrices  B^  are 

=  P(uE)/<°\  (3.37) 

fl(1)  =  -y -V{vE<t>')l (1)  +  P(vB)J<°>  (3.38) 

hx 

and 

B(a)  =  («*(*')’)  ^(3)  -  +  2 v^')/(1)  +  V{v^°\  (3.39) 

The  negative  signs  that  appear  in  the  definitions  of  B^  and  B ^  result  from  the  transposing 
of  /W.  Again,  it  is  noted  that  in  (3.37)  -  (3.39),  the  mx  x  (mx  —  4)  matrices  I^l\l  =  0,1,2 
have  components  <5,-^  as  defined  in  (3.1)  -  (3.3). 

Thus,  the  fourth-order  spatial  problem  (3.6)  can  be  solved  in  a  variety  of  ways  using 
the  Sinc-Galerkin  method.  For  standard  forward  problems,  the  system  (3.16)  is  often  the 
most  convenient  to  formulate  and  solve.  If  the  forward  solver  is  part  of  a  numerical  routine 
for  solving  the  parameter  recovery  problem,  then  (3.28)  is  more  useful  since  El  is  replaced 
by  its  finite  dimensional  approximation.  Both  approaches  yield  solutions  um.  which  are 
exponentially  convergent  approximations  to  the  solution  u  of  (3.6). 

S.2.  The  Spatial  Problem:  Cantilever  Boundary  Conditions 

A  second  set  of  fourth-order  boundary  conditions  arise  when  modeling  beams  that  are 
fixed  at  one  end  and  free  at  the  other.  To  extend  the  Sinc-Galerkin  method  to  problems 
with  these  cantilever  boundary  conditions,  consider  the  ODE 

Lu(x)  =  (EI(x)u"(x))"  =  /(x),  0  <  x  <  1 

u(0)  =  a,  [EIu"){  1)  =  7  (3.40) 

u'(0)  =  J3,  {EIu")'{\)  =  6. 

A  Sinc-Galerkin  method  to  approximate  the  solution  of  (3.40)  can  be  developed  as  follows. 
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Define  the  set  of  basis  functions  by 

B_Mm_A(x),  i  =  -Mx  -  4 
i  =  -Mx  -  3 
B-mw- 2(1),  i  =  —Mx  -  2 
B~Mm-i(x),  i  =  -Mx  —  1 

(i(x)  =  '  v(x)Si(x),  —Mx  <i  <  Nx  (3.41) 

bn.+ i(x),  i  =  Nx  -f  1 

bn.+ i(x),  i  =  Nx  + 2 

BNm+3(,x)>  i  =  Nx  -)-  3 
Bn.+*(x),  i  =  Nz  +  4. 

Here  5,( x)  =  5(i,  A.*)  o  <f>(x)  and  the  basis  weight  v(x)  is  taken  to  be 

v(z)  =  [x(l  -  i)]3.  (3.42) 

The  boundary  basis  functions  are 

H-m.-i(x)  =  (1  -  x)4[20x3  +  10x3  +  4x  +  1], 

bn.+ i(x)  =  x4[20(l  -  x)3  +  10(1  -  x)3  +  4(1  -  x)  +  1], 

B-m.-2(x)  =  x(l  -  x)4[10x3  +  4x  +  1], 
bn,+i{x)  =  -x4(l  -  x)(10(l  -  x)3  +  4(1  -  x)  +  1], 

B-m.~ a(x)  =  x3(l  -  x)4  [2x  +  i  , 

BNm+3{*)  =  X4(l  -  x)3  ^2(1  -  X)  +  ^  , 

B_Mm_, ,(x)  =  ix3(l  -  x)4 

and 

bn.+*(x)  =  -^x4(l  -  x)3. 
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A  brief  note  concerning  the  choice  of  basis  is  in  order  at  this  point.  First,  since 
^y[S(t,  hx)  o  <f>(x)]  ,  l  =  1, 2,  *  •  * ,  are  undefined  at  x  =  0  and  x  —  1,  some  basis  modifi¬ 
cations  must  be  made  when  solving  problems  with  nonzero  boundary  conditions  (see  also 
the  definition  of  tpk  in  (3.23)).  It  is  tempting  to  use  fewer  algebraic  boundary  basis  functions 
and  the  basis  weight 

u(x)  =  x(l  —  x)3, 

but  in  many  problems  this  results  in  nonzero  boundary  terms  when  integrating  by  parts.  By 
using  the  basis  weight  v(x)  =  [x(l  —  x)]3  and  a  full  complement  of  algebraic  terms,  this  pitfall 
can  be  avoided.  Furthermore,  the  basis  {£}  as  defined  in  (3.41)  can  be  used  for  problems 
with  free  boundary  conditions,  thus  providing  consistency  to  the  method. 

The  approximate  solution  is  then  defined  to  be 

um.(x)  =  aC-M.-i(z)  +  /?(-M.-a(*)  +  76v.-t-3(z)  +  $(n.+*(x) 

+U_m.-3C-M.-3(i)  +  ti_M.-4C-Af.-4(®)  +  WAf.  +  lCw,+l(*) 

(3.43) 

UNm-u(Nm+2(x) 

N. 

+  Yj  Ui(i(X) 

i—-Mx 

where  a,/7, 7,  and  6  are  known  and  the  coefficients  {u,}  are  unknown.  The  quantities 

I-eJoF 

and 

j-  It  BI%1) 

~  EI(\)  [£/(1)]>t 

are  well-defined  since  EI(x)  is  assumed  positive  on  [0, 1].  Note  that  with  the  definition  (3.41) 
for  the  basis  {C«})um.(*)  satisfies  the  boundary  conditions;  that  is, 


<.(»)  =  ?,  (£/<J'(l)  =  *• 

The  m,  4-  4  unknown  coefficients  in  (3.43)  are  determined  by  orthogonalizing  the  residual 
with  respect  to  the  set  of  sine  functions  This  Petrov-Galerkin  approach  is  in 
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contrast  to  those  Galerkin  methods  in  which  the  residual  is  orthognalized  with  respect  to 
the  basis  and  is  done  to  take  advantage  of  the  exponential  accuracy  of  point  evaluation  in 
the  quadratures.  This  yields  the  discrete  system 

{Lumm-f,S,)  =  0  (3.44) 

for  p  =  —  Mx  -2,---,Nx  +  2.  The  inner  product  (•,  ■)  is  that  defined  in  (3.8)  with  the  weight 
in  this  case  taken  to  be 

w(x )  =  1.  (3.45) 

The  difference  between  the  weight  function  in  (3.45)  and  that  given  in  (3.9)  is  due  to  the 
presence  of  the  basis  weight  v(x)  in  the  definition  of  the  basis  {(,}.  If  the  definition 

uh(x)  =  52  “»Ci(x) 

i=-M, 

is  made,  then  (3.41)  and  (3.44)  can  be  combined  to  yield 

3)  Sp)u_Mm_3  +  (£5-m.-4>  •5p)u_m.-4  +  (LBNm+i,  Sp)uNm+l 
+(LBNm+ 3,  Sp)utfm+3  +  ( Lu/ „  Sp )  (3.46) 

=  (7,  sr) 

for  p  =  —Ms  —  2,  •  •  • ,  Nx  +  2  and 

J(x)  =  f(x)  ~  aLB_Mm_x{ i)  -  ^L5_w._2( x)  -  iLBNm+3{x)  -  6LBNm+ «(x). 

Integration  by  parts  is  applied  to  ( Lu Sp )  thus  yielding  the  integral 

f  Uh(x)[EI(x)(Sp(x)w(x)y']"dx  (3.47) 

Jo 

(compare  to  (3.10)).  The  weight  choice  (3.45)  is  sufficient  for  guaranteeing  that  the  boundary 
term  (3.11)  vanishes  with  u  replaced  by  «>,.  As  in  Section  3.1  there  are  two  approaches  here. 
In  the  traditional  approach  EI(x )  is  differentiated  directly  and  in  the  alternative  approach 
El  is  simply  replaced  by  the  finite  dimensional  approximate  EImm ■  The  remaining  inner 
products  in  (3.46)  are  evaluated  directly  via  Theorem  2.1. 
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Proceeding  with  the  first  approach,  as  indicated  by  (3.12),  the  choice  of  weight  w  directly 
affects  the  decay  conditions  dictated  by  (2.9)  of  Theorem  2.2.  For  the  weight  v>(x)  =  1  the 
condition 

\EI(x)U(x)\  <  xa+3(l  -  xf +3  (3.48) 

guarantees  sufficient  decay  so  that  the  asymptotic  errors  resulting  from  the  quadrature  can 
be  balanced  by  choosing  hx  and  Nx  as  specified  in  (3.14)  and  (3.15).  The  term  U{x)  denotes 
that  part  of  the  true  solution  which  is  approximated  by  and  is  given  by  the  formal  change 
of  variables 


K{ x)  =  u(x)  -  aB_Mm_i(x)  -  jBNm+3(x) -6BNm+4(x) 

-u"(0)B_M._3(s)  -  u"'(0)B_m._<(x)  (3-49) 

-u(l)Bw.+i(a:)  -  u'(l)BNm+1(x). 


The  discrete  system  for  (3.40)  can  then  be  formulated  as  follows.  Let 
ajv.+i,  aNm+ 2  and  /  denote  the  (m,  +  4)  X  1  vectors  containing  the  product  of  ^  and 
the  approximations  to  the  inner  products  (LB~Mm_3,  Sp)}  (LB_Mm_A,  Sp),  (LB_Nm+i,  Sp), 
(LBffm+2,Sp),  and  (7,  Sp),  respectively.  Hence  the  p-th  entries  of  the  respective  vectors  are 

(*)’  =  w 


(•  =  -M,  -  3,  - Mx  -  4,  AT,  +  1,  N,  +  2)  and 


(!) p  =  7(xp) 


W(Xp) 

<f>'{xpy 


Furthermore,  let  Am  denote  the  ( mx  4-  4)  x  mx  matrix  which  results  from  the  expansion  of 
the  inner  product  ( Luh ,  Sp).  For  a/(x),  t  =  0, 1, 2, 3, 4  as  defined  in  (3.18)  -  (3.22),  the  matrix 
Am  is  given  by 


Am  =  ^I{A)V(vai)  -f  ^I^V(va3)  -f  ^IwV(va2) 

XXX 

+  ^-/(1)P(va!)  +  I^V{vao). 

Here  J<V  =  0,1, 2, 3, 4  are  (mx  +  4)  x  mx  matrices  whose  pt-th  entry  is  given  by  6$  from 
(3.1)  -  (3.5),  and  v  is  defined  in  (3.42). 
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The  discrete  system  for  the  determination  of  the  unknown  coefficients  {u^}  is  given  by 

AXU  =  f 

where  the  (mx  +  4)  x  (mx  +  4)  matrix  Ax  is  defined  to  be 

Ax  =  •  o_a/._3  :  Am  :  a/v.+i  :  oat.+2]-  (3.50) 

Here  tZ  is  defined  to  be  the  (mx  +  4)  x  1  vector 

w  =  +  1,«N.  +  2]T  (3.51) 

containing  the  unknowns. 

It  is  noted  that  the  matrices  Ax  as  defined  in  (3.17)  and  (3.50)  differ  only  in  the  presence 
of  v  in  the  diagonal  multipliers  and  the  addition  of  border  vectors.  Hence  the  method  is  easily 
adapted  when  the  boundary  conditions  are  changed.  Moreover,  the  exponential  convergence 
rate  is  maintained,  thus  preserving  the  accuracy  of  the  method. 

With  parameter  recover  in  mind,  it  is  again  worthwhile  to  use  the  alternative  approach 
to  develop  the  discrete  system  which  arises  when  EI(x)  is  replaced  by  the  finite  dimensional 
term  EImm  as  defined  in  (3.25).  To  simplify  notation  in  the  discussion  which  follows,  let 

N. 

EIh{x)  =  J2  ckVE(x)Sk(x) 

k=-Mm 

and 

EIc(x)  =  C_Af._l5_M.-l(®)  +  C-Af.-2&-M.-2(z)  +  Ctf.  +  l &JV.+1  (®)  +  CNm+7^N.  +  j(x) 

so  that  EImm  =  Elh  +  EIC.  Note  that  EIk  and  EIC  simply  denote  the  sine  and  algebraic 
components  of  the  approximate  parameter. 

Consider  now  the  inner  products  found  in  the  system  (3.46).  The  expansion  of  (Luk,  Sp) 
proceeds  exactly  as  before  with  EI(x)  simply  being  replaced  by  EImm(x)  in  the  integral 
(3.47).  In  the  boundary  inner  products,  ( LB; ,  Sp),i  =  —  Mx  —  4,  ■  •  • ,  —Mx  —  \  and  t  =  Nx  +  1, 
•  ■  •  i  TV*  +  4,  expansion  and  integration  by  parts  yields 

(LBitSp)  =  ||o1EIh(a:)H;'(x)(5pu;y'(x)^  +  BT| 

+  l\EIcB'')"{x)Sp{x)w{x)dx. 

Jo 
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The  weight  w(x)  =  1  (see  (3.45))  is  sufficient  for  guaranteeing  that 

Bt  =  {( EIhB'')'(Spw )  -  EIhBl'(Spw)'}(x)\l  =  0. 

The  resulting  integrals  are  evaluated  via  Theorem  2.2  or,  when  possible,  Theorem  2.1. 

For  the  weight  w(x)  =  1  ,  the  decay  condition  is 

| El{x)U(x)\  <  xa+3(l  -  xf +3 

where  U  and  E 1  are  defined  in  (3.49)  and  (3.27)  (compare  to  (3.48)).  Again,  the  asymptotic 
errors  are  balanced  by  choosing  hx  and  Nx  as  specified  in  (3.14)  and  (3.15). 

The  matrix  system  corresponding  to  (3.40)  may  be  formulated  as  follows.  Let  $(’),  'jd*), 
and  p$(t)  be  defined  as  they  were  in  (3.31),  (3.32),  and  (3.30),  respectively  (with  l  =  0, 1,2 
and  j  =  2,3,4).  Note  that  in  the  definitions  now,  the  index  ranges  are  —  Mx  <  i  <  Nmt 
—  Mx  —  2  <  p  <  Nx  +  2,  and  —Mx  —  2  <  k  <  Nx  +  2,  and  the  (mx  +4)  x  1  coefficient  vector 
is  now  c  =  •  •  • ,  Cffm+2]T-  Hence  $09,  >!>(*)  and  p$(t)  have  the  sizes  (mx  +  4)  x  mx, 

mx  x  (mx  4-  4)  and  mx  x  1,  respectively. 

Furthermore,  let  4"  denote  the  (mx  -f  4)  x  mx  matrix  which  is  defined  componentwise  by 

l*'V  =  ) 

and  let  c  =  [c_m.,  •  •  • ,  c^m]T .  Finally,  for  i  =  — Mx-4,  •  •  • ,  -Mx-1  and  i  =  Nx  +  l,  ■  •  • ,  Nx-f  4, 
let  cLi  denote  the  (mx  +  4)  x  1  vectors 

3,  =  *"V{vEB[')ci-V  (j'(EIcBl')")  f  . 

Here  1  is  simply  the  mx  x  1  vector  of  ones. 

With  the  unknown  vector  u  given  by  (3.51),  the  discrete  system  can  be  written  as 

Axu  =  f 

where 

Ax  =  [a_W(>_4  :  3  :  Am  :  aNa+i  :  3Wb+2]  (3.52) 

and 

/=  £*  *  ~  "  /?3_m._j  -  70W+3  -  SaN+i. 
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The  (rrij.  +  4)  x  mx  submatrix  Am  is  given  by 

Am  =  [^V(vpMi))  +  2  $^D(up*(o)  +  $C4>P(vp*(0,)]. 

It  should  be  noted  that  the  coefficient  matrix  Ax  in  (3.52)  differs  from  that  arising  in 
the  fixed  boundary  problem,  (3.29),  only  in  the  presence  of  v  in  the  diagonal  multipliers  and 
the  addition  of  border  vectors.  This  makes  the  method  easily  adaptable  when  changing  the 
boundary  conditions.  Furthermore,  the  matrices  4>",  and  tyW  can  be  expanded  in  terms 
of  fundamental  matrices  in  a  manner  similar  to  that  in  (3.33)  -  (3.36),  thus  simplifying  the 
implementation  of  the  method.  Finally,  the  exponential  convergence  rate  of  the  method  is 
maintained,  thus  preserving  the  method’s  accuracy. 

With  the  techniques  from  this  section,  the  implementation  of  the  Sinc-Galerkin  method 
for  problems  with  free  and  simple  boundary  conditions  can  be  accomplished  in  a  manner  that 
is  completely  analogous  to  that  used  for  cantilever  boundary  conditions.  Further  details  and 
examples  of  the  Sinc-Galerkin  method  for  problems  with  free  and  simple  boundary  conditions 
can  be  found  in  [8]. 

8.8.  The  Temporal  Problem 

The  last  ODE  to  be  considered  is  the  initial  value  problem 


Pu(t)  =  u(t)  =  / (t),  0  <  t  <  oo 

u(0)  =  tt(0)  =  0. 


(3.53) 


A  Sinc-Galerkin  method  to  approximate  the  solution  of  (3.53)  can  be  developed  in  a  manner 
similar  to  that  of  the  preceding  boundary  value  problems.  Define  the  set  of  basis  functions 

{$;>£-«,  by 

s;(t)  =  s(j,h,)oT(t) 

where  T  :  Dw  — ►  D$  is  given  in  (2.2).  The  approximate  solution  um,(f)  is  then  defined  by 


um,(t)  =  5Z  mt  =  Mt  +  Nt  +  1. 

}--Mt 


(3.54) 


The  mt  unknown  coefficients  {u^}  in  (3.54)  are  determined  by  orthogonalizing  the  residual 
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with  respect  to  the  set  of  sine  functions  which  leads  to  the  analysis  of 

U,s;)  =  (Pu,s;) 

=  («,s;) 

for  q  =  —Mt,‘  •  ■ ,  Nt.  The  weighted  inner  product  for  (3.55)  is  defined  to  be 

(F,G)  =  H  F(t)G(t)w‘(t)dt, 

JO 

and  the  weight  is  taken  to  be 


(3.55) 


w*(0  =  ~n 


1 


JW) 

for  reasons  that  are  discussed  in  [4].  As  before,  integration  by  parts  is  used  to  transfer  the 
differentiation  of  u  onto  S*w*.  To  guarantee  that  boundary  terms  vanish,  it  is  assumed  that 

lim  -  t-  =  0. 

*-*o  ln(t ) 

The  resulting  integrals  are  then  evaluated  via  the  quadrature  rules  of  Section  2.  With  respect 
to  (2.4),  the  condition 

n —  ,  I  <€(0,l) 

KOVmOI  ^  ft  | 

[  t's,  t  G  [l,oo) 

guarantees  the  boundedness  necessary  to  truncated  the  infinite  quadrature  rule.  With  7  and 
6  specified  and  Mt  chosen,  the  parameter  selections 


ht  = 


I  ltd 
7  Mt 


and 


Nt  =  |  lMt  +  lj  (3.56) 

balance  the  asymptotic  errors  in  (2.5)  to  at  least  0(e~(1rd'rM,^ ). 

In  many  time-dependent  PDE’s,  it  is  reasonable  to  assume  that  the  solution  decays 
exponentially  at  infinity;  that  is,  the  solution  satisfies 

l»W\/t(0l  <  * 
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t\  t  e  (0,1) 

e~St,  t  G  [1,  00) 


or,  more  stringently, 


(3.57) 


|u(f)|  <  KP+h-“. 

With  this  supposition,  Lund  [12]  shows  that  the  condition  (3.56)  can  be  replaced  by 

Nt  =  |^-ln  +  l]  .  (3.58) 

The  selection  Ni  in  (3.58)  significantly  reduces  the  size  of  the  discrete  system  with  no  loss 
of  accuracy. 

The  discrete  system  for  (3.53)  can  then  be  formulated  as  follows.  Let  I^l\l  =  0,2  denote 
the  mt  x  mt  matrices  whose  qj- th  entry  is  6$  from  (3.1)  and  (3.3),  and  let  T>(r])  again  be 
the  diagonal  matrix  with  entries  t?(t_ at,),  •  •  • ,  T]{^Nt)-  With  the  usual  definitions  for  U  and  f 
and  the  identity 

(f(())-*{(T(())-i}"  =  -i, 

the  system  for  the  determination  of  the  unknown  coefficients  {ifj}  is  given  by 

Atu  =  Z’((T)~*)/  (3.59) 

where 

A  =  i/(,)  -  \l  »((*)*).  (3.60) 

Further  details  concerning  the  derivation  of  the  system  (3.59)  can  be  found  in  [4]. 

Note  that  nonzero  initial  conditions  can  be  handled  in  a  manner  analogous  to  that  used 
for  nonzero  boundary  conditions  in  the  previous  discussion.  Rational  initial  basis  functions 
are  used  to  incorporate  the  initial  behavior  and  this  known  contribution  is  then  taken  to  the 
right-hand  side  of  the  resulting  discrete  system  (see  [8]).  All  other  analysis  is  identical  to 
that  for  the  problem  (3.53). 

4  Time-Dependent  Problems 

This  section  details  the  Sinc-Galerkin  method  applied  to  fourth-order  time-dependent  PDE’s 
with  fixed  and  cantilever  boundary  conditions.  Since  the  choice  of  basis,  test  functions,  and 
inner  product  are  all  straightforward  extensions  of  those  used  to  solve  the  ODE’s  in  Section  3, 
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the  error  analysis  and  system  formulation  follow  directly  from  previously  discussed  results. 
Once  a  discrete  system  has  been  formulated,  various  options  exist  for  solving  the  associated 
matrix  equation.  Two  such  algorithms  are  outlined  and  their  relative  merits  for  various 
problems  are  discussed. 

4-1-  The  Time-Dependent  Problem:  Fixed  Boundary  Conditions 
Consider  the  time-dependent  problem 

d7u  d 2  (  d2u  \ 

“  dt’C*’*)  +  0?  [EI(x)fcj(xit)]=f(xit)'  0<x<l  t  >  0 

u(0,  t)  =  u(l,  t)  =  0,  t  >  0 

!«».*) -£(>.*)« o.  *>«  <*■« 

du 

u(x,0)  —  —  (x,0)  =  0,  0  <  x  <  1. 

Given  the  basis  {S^}  where 

Sij(xtt)  =  Si(x)s;(t)  =  S{i,  hx)  o  <f>(x)S(j,  ht)  O  T(t), 

the  approximate  solution  is  defined  by  way  of  the  tensor  product  expansion 

N.  Nt 

uySy(*.0»  mx  =  Mx  -f  Nx  +  1,  mt  =  Mt  +  Nt  +  1. 

t  —  —  Mm  | 

The  mx  ■  mt  unknown  coefficients  {u,/}  are  determined  by  orthogonalizing  the  residual  with 
respect  to  the  set  of  sine  functions  This  yields  the  discrete  Galerkin 

system 

-  f,SpS *)  =  0 

for  p  =  —  Mx,  ■  ■  ■ ,  Nx  and  q  =  —  Mf,  •  ■  • ,  Nt.  The  inner  product  (•,  •)  is  taken  to  be 

(F,G)  =  f  f  F(x,t)G(x,t)w(x,t)dxdt  (4.2) 

Jo  Jo 

where 

w(x,t)  =  w(x)u;*(0  =  (^(x))_?(t(f))_^. 

The  quadrature  rules  and  one-dimensional  results  from  Sections  3.1  and  3.3  can  be  used  to 
determine  the  resulting  matrix  system. 
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As  before,  equating  asymptotic  errors  is  fundamental  to  minimizing  system  size.  When 
the  decay  conditions  (3.13)  or  (3.26),  and  (3.57)  are  combined  to  yield 

\El(x)u(x,t)\  <  Kxa+l(\  -  xf+ip+le-Si  (4.3) 


then  the  choices 


\£l(x)u{x,t)\  <  Kxa+I(  1  - 


*■  -  fe 

hi  ~  hX) 


NX=  JgM.+  lJ, 


Mt=  -Mx  + 1  ,  (4.8) 

11.7  JJ 

and 

(4-9) 

for  the  stepsizes  and  summation  limits  balance  the  asymptotic  errors.  If  one  takes  d  — 
then  the  above  choices  yield  an  asymptotic  error  rate  of  order  0(e~"\/'aMm/2). 

Given  Mx,  Nx,  Mt,  Nt  and  h  =  hx  =  ht  as  defined  above,  the  discrete  system  for  (4.1)  is 


AXUCJ  +  CXUA[  =  G. 


(4.10) 


°-mV\7  ■ 


_  f  XV  \  /w*\ 

G  =  V\-)  FV(-r-). 

\<t>')  VT  J 


(4.11) 

(4.12) 


The  mt  x  mt  matrix  At  is  given  by 


*  . 


(4.13) 
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as  shown  in  (3.60).  Furthermore,  mx  x  mt  matrices  U  and  F  are  defined  componentwise  by 


W)ij  =  Uij 

and 

[Fhj  = 

It  should  be  noted  that  the  ordering  of  the  coefficients  in  U  mimics  that  used  in  most 
standard  time-differencing  schemes.  This  is  a  matter  of  convenience  since  the  Sinc-Galerkin 
method  is  not  bound  by  any  specific  ordering  of  the  grid. 

The  structure  of  the  m,  x  mx  matrix  Ax  depends  on  the  scheme  that  is  used  to  discretize 
EI(x).  If  the  parameter  is  fully  differentiated,  then  Ax  is  given  by  (3.17).  If,  on  the  other 
hand,  EI(x )  is  approximated  by  a  linear  combination  of  sine  and  algebraic  basis  functions, 
then  Ax  is  given  by  (3.29). 

Various  methods  exist  for  solving  the  equation  (4.10).  Referred  to  as  a  generalized 
Sylvester  equation  (4.10)  is  algebraically  equivalent  (page  414  of  [15])  to  the  system 

An  =  { Ct  ®AX  +  At®  Cx)  co(U)  =  co(G )  (4.14) 

where  the  tensor  or  Kronecker  product  of  an  m  x  m  matrix  E  with  a  p  x  q  matrix  H  is 
defined  by 

E  ®  H  =  [ eijH}mPXnq ■ 

—4 

The  vector  u  =  co(U )  is  the  concatenation  of  the  mx  x  mt  matrix  U  obtained  by  successively 
“stacking”  the  columns  of  U,  one  upon  another,  to  obtain  an  mxmt  x  1  vector. 

The  system  (4.14)  can  be  solved  directly  via  any  of  the  decomposition  methods  that  are 
available  for  linear  systems.  Although  this  system  is  easily  formulated,  the  fact  that  A  is 
very  large  (mz  •  mt  x  m,  •  mt)  and  not  banded  causes  this  method  to  be  impractical  in  some 
problems.  For  more  general  fourth-order  operators  however,  this  may  be  the  only  method  of 
choice. 

A  second  algorithm  for  solving  (4.10)  depends  on  the  generalized  Schur  decomposition 
(page  396  of  [16]).  As  guaranteed  by  the  results  of  Moler  and  Stewart  [17],  there  exist  unitary 
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matrices  Qi,ZifQ2,  and  Z2  such  that 


Q\AXZX  =  P 
Q\CXZX  =  R 
Q\GtZ2  =  S 
Q\AtZ2  =  T 


where  P,R,S,  and  T  are  upper  triangular.  If  Y  =  Z*UZ2  and  C  =  Q[GQ2,  then  (4.10) 
transforms  to 


PYT *  +  RYS *  =  C. 


By  comparing  the  A:-th  columns,  one  finds  that 


n  n 

P  J2  ik]Vj  +  R  X]  SkiVi  =  °k 

i—k  i-k 


which  yields 

(tkkP  +  skkR)yk  =  cfc  -  P  ikjVj-R  Y2  Sk>yj  (415) 

j=k+ 1  j'=fc+ 1 

(for  convenience,  it  is  assumed  that  all  matrices  are  n  x  n  and  indexed  from  1  to  n).  With 
the  assumption  that  the  matrix  (tkkP  +  skkR )  is  nonsingular,  the  solution  to  (4.15)  is  easily 
found  by  recursively  solving  triangular  systems. 

Although  this  algorithm  does  require  complex  algebra,  it  is  both  robust  and  efficient 
and  requires  no  assumptions  concerning  the  diagonalizability  of  the  component  matrices. 
It  should  be  no  ed  that  a  “real”  version  of  this  algorithm  also  exists  (18).  In  this  latter 
algorithm,  Qi,Zx,Q2,  and  Z2  are  orthogonal  with  P,S  quasi- upper  triangular  and  R,  T 
upper  triangular. 
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4.2.  The  Time-Dependent  Problem:  Cantilever  Boundary  Conditions 

A  generalization  of  the  problem  which  arises  when  modeling  beams  with  cantilever  bound¬ 
ary  conditions  is 

d2u . 


Cu(x,t)  = 


EI(x)-^{x,t)j  =  f{x,t),  0<x<l 


x’t)+b{ 

u(0,  t)  =  a(t),  (M)  =  7(0»  *  >  0 

(o,t)  =  5(i),  k(«5?)<M>-*W. 


dx 


(4.16) 


5u , 


«(*,°)  =  =  °> 


t  >  0 


0  <  x  <  1. 


The  basis  for  this  problem  is  taken  to  be  {£/(®)SJ(<)}  where  Ci(x)  1®  defined  in  (3.41)  and 
Sj(t)  =  S(j ,  ht )  o  T(<).  Here  the  approximate  solution  is  taken  to  be 


u 


mtmi 


(*.0  =  53  53 

t= — Mb  j—  -Aft 
N« 

+  53  <^i(0{u-W.-3jC-W.-3(*)  +  W-M.-4,jC-M.-4(*) 
j=-Mt 

+urt.+i,j(Nm+i(x)  +  u^.+2,jCAr.+a(*)} 

+{a(t)C_M.-i(»)  +  ?C-M.-a(*)  +  7(<Kfr.+3(*)  +  ^(0Cw.+4(*)} 


where 


and 


7(0  = 


E/(l) 


7(0 


*(0  = 


1  EI'{  1) 


-?(0  - 


7(0- 


EI(  1)  v  7  [£/(l)]a 

It  should  be  noted  that  the  approximate  solution  does  satisfy  the  boundary  conditions  in 
(4.16). 

The  (mx  +  4)  •  mt  unknowns  {ttj;}  are  determined  by  orthogonalizing  the  residual  with 
respect  to  the  sine  functions  {Sp{x)Sl(t)}~^M^t2  .  Na+J.  This  yields  the  discrete  system 


( Cum,mt  -  f,SpS *)  =  0, 
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-Mx  -2<p<Nx  +  2 
-Mt  <q<Nt 


where  (•,  •)  is  defined  in  (4.2)  with  w(x,t )  =  (T(J))“i. 

Appropriate  integration  by  parts  and  the  application  of  the  one-dimensional  results  from 
Sections  3.2  and  3.3  yields  the  matrix  equation 

AxUCj  +  CxMUAj  =  G  (4.17) 

where  Ct,Cx  and  At  are  defined  in  (4.11),  (4.12),  and  (4.13)  respectively,  and 

The  (m,  +  4)  x  mt  matrices  U  and  F  are  defined  componentwise  by 


and 

=  7{Xii  tj) 

where 

7(*,0  =  /(*.t)-£(a(0fl.«..,(x))-£(]5(0B.Jl,..,(*)) 
-CW)Br.M)  -  C(S(t)B„.+ ,(*)). 

The  (mx  -f  4)  X  (mx  -f  4)  matrix  Ax  is  given  by  (3.50)  or  (3.52)  depending  on  which  scheme 
is  used  to  discretize  El.  Finally  the  (mx  -f  4)  x  (mx  +  4)  matrix  M  has  the  form 


where  the  mx  x  mx  submatrix  Dm  and  the  ( mx  4-  4)  x  1  vectors  are  giver;  by 


Dm  =  D(v), 


i>L2  =  D(B-M.-*)1) 
i>Ll  =  D(B_Mm- s)f, 
=  D(Bnm+ i)f 
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and 


bR2  =  V(Bfijm+ 2)L 

The  system  (4.17)  can  again  be  solved  via  the  generalized  Schur  algorithm  (4.15)  as  discussed 
in  Section  4.1. 

Before  implementing  the  method,  the  decay  parameters  6,/3, 7  and  6  must  be  determined 
and  summation  limits  chosen.  For  the  spatial  weight  w(x )  =  1,  the  decay  conditions  are 

| El{x)U(xt  f)|  <  xa+3(l  -  xf+3r+ie-st  (4.18) 

or 

\£l(x)U{x,t)\  <  xa+3(l  -  x)*+3r+*e-“  (4.19) 

depending  on  the  manner  in  which  El  was  discretized.  Here  £X(x)  is  defined  in  (3.27)  and 
U(x,t)  denotes  that  part  of  the  true  solution  which  is  approximated  by 

Mx>t)=  E 

»= —  M.  j=-M( 

(see  also  (3.49)).  With  the  decay  parameters  specified  and  M*  chosen,  the  remaining  stepsizes 
and  summation  limits  are  given  by  (4.5)  -  (4.9). 


5  Numerical  Examples 

The  four  examples  reported  in  this  section  were  selected  from  a  large  collection  of  problems 
to  which  the  Sinc-Galerkin  method  was  applied.  The  results  are  representative  of  those 
obtained  on  other  sample  problems.  For  purposes  of  comparison,  contrast,  and  performance 
evaluation,  examples  with  known  solutions  were  chosen.  The  first  and  last  examples  have 
analytic  solutions,  the  second  example  has  an  algebraic  singularity,  and  the  third  example 
contains  a  logarithmic  singularity.  As  will  be  demonstrated  by  the  numerical  results,  the 
boundary  singularities  have  no  adverse  affect  on  the  performance  of  the  method. 
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In  all  examples  d  =  ?.  The  errors  are  reported  on  both  the  set  of  sine  gridpoints 


and  the  set  of  uniform  gridpoints  (lx  =  j$,lt  = 


_ 


U  =  {(2p.  j«)}J=!;‘.‘.'.b8o°  i  2p  =  t*P>  3o  = 


The  errors  on  these  grids  are  reported  as 


\ES  (^x,/lf)||  —  ^i)  ~  wmimt(®i>  tj)\ 

—  M b  jV( 


W^U  lU(2P>  3<l)  um,m«(zp>  5?)l> 

0<?<50 

respectively.  The  dependence  of  both  errors  on  the  number  of  sine  basis  functions  is  indicated 
with  the  superscript  Mx.  It  is  noted  that  if  exponential  convergence  is  realized,  then 

where  Mx  and  Mx  denote  the  lower  limits  for  the  spatial  sums.  In  the  examples  of  this  section, 
Mx  =  2 Mx,  and  the  exponential  convergence  is  verified  by  comparing  (|| E^"^2(£X}  ^f)||)v^  and 

The  error  and  convergence  results  are  tabulated  in  the  form  .aaa  —  7  which  represents 
.aaa  x  lO-7.  All  problems  were  run  with  sixteen  place  accuracy  on  a  Vax  8550. 


Example  5.1. 

d’i  u  d 1  (  d7u  \ 

+  dx*  \  ~  0<x<1  i>0 

tt(0,i)  =  u(l,t)  =  0,  t  >  0 

5(0.0-50.0-0.  <>o 

du 

Ti(x,0)  =  — (cc,0)  =  0,  0<x<l 

The  flexural  rigidity  is  EI(x)  =  14-  sin(7rx)  and  /(x,  t)  is  consistent  with  the  true  solution 
u(x,t)  =  x(l  -  x)8in(47rx)tJe“‘.  The  parameter  is  expanded  via  (3.25)  thus  yielding  the 
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spatial  representation  (3.29).  The  decay  condition  (4.4)  dictates  the  parameter  choices  a  = 
/3  =  7  =  |  and  S  =  1.  The  asymptotic  error  rate  0(e~*^3M’/3)  is  maintained  as  indicated 
by  the  last  column  of  Table  1.  Notice  that  the  choice  of  Nt  given  in  (4.9)  significantly 
reduces  the  size  of  the  matrices  involved  in  (4.10).  The  errors  on  both  the  sine  grid  and 
the  uniform  grid  do  not  differ  dramatically  though  in  the  next  example  the  difference  will 
be  more  noticeable.  Figure  3  shows  the  true  solution  u(x,  t)  while  Figure  4  shows  both  the 
true  and  approximate  solutions  (for  Mx  —  8, 16)  at  the  time  slice  t  =  2.  The  approximate 
solution  for  Mx  =  32  is  buried  in  the  true  solution  on  this  scale. 


Mx 

Nx 

Mt 

N,  Ml 

lK-(^.)ll 

4 

4 

4 

2  .6207  -  0 

.8040  -  0 

8 

8 

8 

4  .6890  -  1 

.8448  -  1 

.7345  -  0 

16 

16 

16 

6  .2910  -  3 

.1105-2 

.3035  -  1 

32 

32 

32 

9  .1906  -4 

.2151-4 

.6587  -  4 

Table  1.  Errors  on  the  Sine  Grid  S  and  the  Uniform  Grid  U  for  Example  5.1. 
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x  axis 


Figure  4.  True  and  Approximate  Solutions  to  Example  5.1  at  Time  t  =  2 
- (Af,  =  8),  -  •  -  (M*  =  16),  - (True). 
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Example  5.8. 

ffiu  q*u 

~Qp(.x'  0  +  fai(x' 0  =  /(*»  0.  0  <  i  <  1  t  >  0 

u(0,  t)  =  u(l,t)  =  0,  t  >  0 

|(°,0=|(M)  =  °,  oo 

«(*»0)  = -gj^C*.0)  =  °.  0  <  x  <  1 

The  function  f(x,t)  i8  consistent  with  the  solution  u(x,t)  =  [x(l  —  x)]7^t6^7e~t  which  has 
algebraic  singularities  at  x  =  0,  x  =  1  and  t  =  0.  The  spatial  discretization  is  taken  to  be 
(3.17)  with  the  decay  parameters  a  =  /?  =  7  =  2, 6=1  dictated  by  (4.3).  As  indicated  by 
Table  2,  the  asymptotic  error  rate  0(e-,rv^7’)  is  achieved  in  spite  of  the  boundary  singu¬ 
larities.  The  increased  accuracy  of  the  method  for  this  problem  as  compared  to  Example 
5.1  is  due  to  the  larger  values  of  a,/3  and  7.  Here  the  error  on  the  sine  grid  is  substantially 
smaller  than  that  on  the  uniform  grid.  This  emphasizes  that  one  should  use  caution  when 
assessing  performance  of  a  method  based  on  only  the  errors  at  the  gridpoints  of  the  method. 
The  true  solution  is  plotted  in  Figure  5  while  time  slices  (at  time  t  —  2)  of  both  the  true 
and  approximate  solutions  (for  Mx  =  4,8)  are  plotted  in  Figure  6. 


Mm 

N. 

Mi 

Nt 

4 

4 

4 

3 

.2040  -  4 

.3184-3 

8 

8 

8 

4 

.9361  -  5 

.3618  -  4 

.1134-4 

16 

16 

16 

7 

.1372  -6 

.1788  -5 

.5233  -  6 

32 

32 

32 

11 

.1742  -  9 

.1338  -  7 

.7441  -  8 

Table  2.  Errors  on  the  Sine  Grid  S  and  the  Uniform  Grid  U  for  Example  5.2. 
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Figure  5.  True  Solution  4,0  Example  5.2. 


*10' 


Figure  6.  True  and  Approximate  Solutions  to  Example  5.2  at  Time  t  =  2 
- (Mx  =  4),  -  •  -  (Mx  =  8),  - (True). 
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Example  5.S. 

d*u  d*u 

-gp(x>  *)  +  0  =  /(* *  0.  0  <  X  <  1  t  >  0 

u(0,t)  =  0,  ~(l,t)  =  2iae-'J  i  >  0 

du  d3u 

^(o,<)  =  o,  ^(1,0  =  0,  (>o 

du 

u(x,0)  =  — (x,0)  =  0,  0  <  X  <  1. 

C/l 

The  true  solution  u(x,t)  =  [(x/n(x))4  +  xJ]tJe_<  dictates  the  forcing  function  /(x,<).  The 
decay  condition  (4.18)  with  EI(x )  =  1  yields  the  parameter  choices  a  =  (3  =  1,7  =  |, 
and  6  —  1  which  in  turn  implies  the  asymptotic  error  rate  0{e~r^M‘^2).  As  indicated  by 
the  results  in  Table  3  this  rate  is  achieved  in  spite  of  the  logarithmic  singularity  at  x  =  0. 
The  convergence  of  the  method  is  even  accelerated  which  can  be  seen  by  in  the  last  column 
of  Table  3.  The  mesh  plot  in  Figure  7  shows  the  distinctive  behavior  that  the  solution  can 
exhibit  when  cantilever  boundary  conditions  are  in  force.  The  “oscillation”  at  the  right-hand 
end  is  tracked  accurately  by  this  method.  The  time  slice  at  t  —  2  shown  in  Figure  8  shows 
the  approximate  ( Mx  =  4,8, 16)  as  well  as  the  true  solution. 


Mx 

Nx 

Mt 

Nt 

l|E?-(*>,A)ll 

ii#*(4,«n 

(l|£"-/1(4,<,)||)'/5 

4 

4 

3 

2 

.4227  - 1 

.9083  -  1 

8 

8 

6 

3 

.8034  -  2 

.1999  -  1 

.3363  -  1 

16 

16 

11 

4 

.8799  -  3 

.1838  -2 

.3953  -  2 

32 

32 

22 

7 

.3947  -  4 

.8741  -4 

.1353  -2 

Table  3.  Errors  on  the  Sine  Grid  S  and  the  Uniform  Grid  U  for  Example  5.3. 
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I 


Figure  7.  True  Solution  to  Example  5.3. 


■~>ximate  Solutions  to  Example  5 


Figure  8.  True  and 


UllKl] 


Example  5.^. 


d7u 


d 2 


d2u 


(*.0  +  ^i(^(*fc(*.0|=/(*.0.  o  <  i  <  i  t>o 


dt> 


'dx* 


u(o,  t)  =  o,  (1.0  =  0,  *>o 


A 

u(x,0)  =  — (®,0)  =  0,  0  <  x  <  1. 

(/t 


3j2--‘  t  >  0 


This  example  illustrates  a  problem  where  the  spatial  discretization  (3.52)  is  useful.  Here 
the  flexural  rigidity  is  EI(x)  =  1  +  sin(Trx)  and  the  forcing  function  f(x,t )  is  consistent 
with  the  true  solution  u(x,t)  =  [sin3(7ri)  +  x]t2e~l.  The  parameter  EI(x)  is  expanded  via 
(3.25)  thus  yielding  the  spatial  representation  (3.52).  The  parameter  choices  a  =  f3  =  1, 
7=|  and  5=1  follow  from  (4.19).  As  demonstrated  by  the  last  column  of  Table  4,  the 
asymptotic  error  rate  0(e_irv^WJ)  |B  achieved  for  larger  values  of  Mx ,  Nx,  Mt,  and  Nt.  On 
this  example  the  errors  on  the  sine  grid  and  those  on  the  uniform  grid  are  nearly  the  same. 
The  smaller  parameters  a  and  indicate  why  the  errors  here  are  larger  than  in  the  previous 
three  examples.  A  mesh  plot  of  the  true  solution  is  shown  in  Figure  9  while  a  time  slice 
( t  =  2)  of  both  the  true  and  approximate  solutions  (Mx  =  8, 16)  are  plotted  in  Figure  10. 


Mx 

Nx 

Me 

Nt 

\\Ep{hm,ht)\\ 

ii  *#•(<.,  <.)« 

8 

8 

6 

3 

.6889  -  1 

.6958  -  1 

16 

16 

11 

4 

.2958  -  1 

.3165-1 

.2307  -  1 

32 

32 

22 

7 

.4346  -  3 

.4885  -  3 

.7572  -  2 

Table  4.  Errors  on  the  Sine  Grid  S  and  the  Uniform  Grid  U  for  Example  5.4. 
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0  0  1  0  2  0  3  0  4  0  5  0  6  0  7  0  8  0,9 

%  a*i» 


Figure  10.  True  and  Approximate  Solutions  to  Example  5.4  at  Time  t 
- (Mx  =8,  )  -  •  -  (A/*  =  16),  - (True). 
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6  Conclusions 


A  fully  Sinc-Galerkin  method  in  both  space  and  time  is  presented  for  fourth-order  time- 
dependent  problems  with  fixed  and  cantilever  boundary  conditions.  The  sine  basis  properties 
which  facilitate  the  simple  assembly  of  the  discrete  system  are  discussed  in  Section  2.  In 
Section  3,  the  sine  discretizations  for  the  second-order  temporal  problem  and  the  fourth-order 
spatial  problems  are  presented.  Alternate  formulations  for  the  variable  parameter  fourth- 
order  problems  are  given  which  prove  to  be  especially  useful  when  applying  the  forward 
techniques  of  this  paper  to  parameter  recovery  problems.  The  ODE  results  are  then  combined 
in  Section  4  to  form  the  discrete  systems  corresponding  to  the  time-dependent  problems  of 
interest.  Computational  issues  are  discussed  and  a  robust  and  efficient  algorithm  for  solving 
the  resulting  matrix  system  is  outlined.  Numerical  examples  which  highlight  the  method  are 
given  in  Section  5.  As  demonstrated  by  the  numerical  results,  the  exponential  convergence 
rate  of  the  method  is  maintained  for  problems  with  both  analytic  and  singular  solutions  as 
well  as  fixed  and  cantilever  boundary  conditions. 
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